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Abstract 

The Lennard-Jones (LJ) Potential Energy Problem is to construct the most 
stable form of N atoms of a molecule with the minimal LJ potential energy. This 
problem has a simple mathematical form 



N N 

minimize 



/(a;) ^ (^^^j subject to a; e M", 

where = (x3,_2 - a;3j_2)^ + (xs^^i - xsj-i)^ + {x^i ~ x^j)'^, (a;3i_2, xa^-i, xa^) 
is the coordinates of atom i in R^, i,j — 1,2, .. . ,N{> 2 integer), and n = 3N; 
however it is a challenging and difficult problem for many optimization methods 
when N is larger. In this paper, a brief review and a bibliography of important 
computational algorithms on minimizing the LJ potential energy are introduced 
in Sections 1 and 2. Section 3 of this paper illuminates many beautiful graphs 
(gotten by the author nearly 10 years ago) for the three dimensional structures of 
molecules with minimal LJ potential. 

1 Introduction 

the Lennard-Jones potential energy problem is to construct the most stable form of N 
atoms of some material with the minimal energy structure. Its form in mathematics is 
very simple: 



minimize f{x) = 4 > > \ ~a T subject to x G M", (1) 




where Tij = (x3i_2 - 3:3^-2)^ + (a^si-i - x^j-i f + {xsi - x^j)'^, {xsi_2, x^i-i, xsi) e 
is the coordinates of atom i, i,j = 1,2, .. . ,N(> 2 integer), and n = 3N. However, 
it interests many researchers in the field of biology, physics, chemistry, mathematical 
optimization, computer science, and materials science. The reason lies in the noncon- 
vexity of the objective function and the huge number of local minima, which is growing 
exponentially with N. [22] tells us the number of distinct local minima of A^-atoms 
Lennard-Jones problem is about 0{e^ ). For example, When =6, 7, 8, 9, 10, 11, 
12, 13 the numbers of distinct local minima of the problem are 2, 4, 8, 18, 57, 145, 
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366, 988 respectively ([26j)- From the point of view of numerical optimization meth- 
ods, this problem is an excellent test problem for local or global optimization methods. 
However, the Lennard- Jones Problem is very difficult and challenging to many opti- 
mization methods even when N is not very large. By now, the best objective function 
vlues for 2 < iV < 309 [Ml [Sni EI] are well known. At first, around 1972 Hoare and 
Pal [231 [Ml [25] gave the best values for atoms 5~16, 18~21, 25, 26, 29, 55. Then, 

in 1987, Northby ([35j) made a landmark yielding most of the lowest values in the 

range 13 < < 150; and in 1999 [39] presented the best results on 151 < A < 309 
with the exception of A =192, 201 [l3], A^4=200, 300 [36], A =185 [3l], A^ =186, 
187 [19j . At present, about thirty percent of Northby's results have been improved 
[46l[l2l[l3l[32l[l6l[l5l[44l[9l[IIl[l5]. The improvements are described as follows. In 
1994 Xue ([16]) got the best values for A^ =65,66,134. Doye et al [l3j) in 1995 
gave the best values on A=38,75~77,102~104; and Leary in 1997 ([32]) presented best 
values for A =69,78,88,107,113,115. Wolf ([l5]) got the best values on A =82, 84, 86, 
88, 92, 93, 94, 95, 96, 99, 100 in 1998. For the left so-called magic numbers sequence: 
17,23,24,72,88, Freeman, Farges, Wille, Coleman, Deaven ( [Ml [IH [Ml [9l [H] ) got the 
best results, respectively. 

The outline of this paper is as follows. In section 2, we review the methods in the 
most important references above-mentioned. Section 3 contains a description of our 
techniques and methods to tackle problem ([TJ; and results of numerical experiments 
for the problem are given in this section. The structure of our optimal solutions is 
illustrated, and can be compared with those of others, at the end of this section. 
Section 4 concludes the paper. 

2 Approaches to the Lennard-Jones Problem 

Some good summaries on the methods solving the Lennard-Jones problem in fact can 
be found in [35l [IH [361 US [32] . In this section we review the methods in the most 
important references [23l [Ml [SSI [35l [Ml llll [46l [111 [11 [121 [ill [^ 

Hoare and Pal's work ([231 IMl [25]) may be the early most successful results on 
Lennard-Jones problem. The idea is using build-up technique to construct the initial 
solutions which are expected to represent low energy states, and using those initial 
solutions as starting points for a local search method to relax to the optimal solution 
([25j). The starting seed is the regular unit tetrahedron with atoms at the vertices, the 
obvious global optimal solution for A = 4. Beginning with this tetrahedron, Hoare and 
Pal added one atom at a time to construct a sequence of polytetrahedral structures 
and at last got good results up to A = 66. For example, for A = 5 its globally 
optimal trigonal bi-pyramid (bi-tetrahedron) structure is gotten by adding an atom at 
the tetrahedral capping position over a triangular face; following the bi-tetrahedron 
structure, the optimal structure of A = 6 is tri-tetrahedron (another known optimal 
structure for A = 6 is octahedron (using tetrahedral capping over triangular faces and 
half-octahedral capping over square faces), which is not a polytetrahedron); for N = 7 
its best structure constructed is the pentagonal bi-pyramid, a structure with a five-fold 
axis of symmetry. Many computer science data structure procedures such as greedy 
forward growth operator and reverse greedy operator can make the build-up technique 



2 



work well. [23j describes the application of methods of studying noncrystalline clusters 
to the study of "spherical" face centred cubic (fee) microcrystallites. [23] describes the 
chief geometrical features of the clustering of small numbers of interacting particles. 

The data structure of Northby (^35j) in finding the good starting solution is the 
lattice based structure. The lattice structures consist of an icosahedral core and partic- 
ular combinations of surface lattice points. |34j first constructed a class of icosahedral 
packings by adding successively larger icosahedral shells in layers around a core central 
atom; this icosahedral lattice can be described as 20 slightly flattened tetrahedrally 
shaped fee units with 12 vertices on a sphere centered at the core atom. Atoms within 
each triangular face are placed in staggered rows in a two dimensional hexagonal close- 
packed arrangement. Each atom in the interior of a face in a given shell is a tetrahedral 
capping position relative to three atoms in the underlying shell. Northby relaxed the 
structure of |34j to get his IC and FC multilayer icosahedral lattice structures. The 
IC lattice can be referred to the FORTRAN code in [46j; it consists of all those sites 
which will comprise the outer shell of the next complete Mackay ([34]) icosahedron. FC 
lattice is a slight modification of IC lattice in that its outer shell maintains icosahedral 
symmetry and consists of points at the icosahedral vertices and the stacking fault posi- 
tions of the outer IC shell. Basing on the IC and FC lattices, Northy gave his algorithm 
first finding a set of lattice local minimizers and then relaxing those lattice minimizers 
by performing continuous minimization starting with those lattice minimizers. The 
algorithm was summarized as Algorithm 1 and Algorithm 2 of |46j . 

The great majority of the best known solutions of Northy are icosahedral in char- 
acter. The hybridization of global search and local search methods, usually, is more 
effective to solve the large scale problem than the global search method or local search 
method working alone. Catching those two ideas, [39] combined a genetic algorithm 
with a stochastic search procedure on icosahedrally derived lattices. The structures of 
the optimal solutions gotten in [39] are either icosahedral or decahedral in character. 
The best results of [45j for iV = 82,84,86,88,92,93,94,95,96,99,100 were gotten by 
using a genetic algorithm alone. Deaven et al ([H]) also using the genetic algorithm 
got the optimal value known for the magic number N = 88. 

The successful works to improve Northby's results ([35]) were mainly done by Xue 
([16]), Leary ([32]), and Doye et al ([13 [l3]). 

Xue ([46]) introduced a modified version of the Northby algorithm. He showed 
that in some cases the relaxation of the outer shell lattice local minimizer with a worse 
potential function value may lead to a local minimizer with a better value. In Northby's 
algorithm the lattice search part is a discrete optimization local search procedure, 
which makes a lattice move to its neighboring lattice with 0{N3) time complexity. 
In [17] Xue introduced a simple storage data structure to reduce the time complexity 
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to 0{Ns) per move; and then used a two-level simulated annealing algorithm within 
the supercomputer CM-5 to be able to solve fastly the Lennard-Jones problem with 
sizes as large as 100,000 atoms. In [l6] by employing AVL trees (^|) data structure 
Xue furthermore reduced the time complexity to 0(log N) if NN (nearest neighbor) 
potential function is used. He ([M]) relaxed every lattice local minimizer found instead 
of relaxing only those lattice local minimizers with best known potential function value 
by a powerful Truncated Newton local search method, and at last got the best results 
known for N = 65, 66, 134, 200, 300. 
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Leary ([32]) gave a successful Big Bang Algorithm in 1997 for getting the best values 
known of iV = 69,78,88,107,113,115. In [32] the FCC lattice structure is discussed 
and its connections are made with the macrocluster problem. It is also concluded in 
[32| that almost all known exceptions to global optimality of the well-known Northby 
multilayer icosahedral conformations for microclusters are shown to be minor variants 
of that geometry. The Big Bang Algorithm contains 3 steps: Step 1 is an initial solution 
generating procedure which randomly generates each coordinate of the initial solution 
with the independently normal distribution; Step 2 is to generate the new neighborhood 
solution by discrete-typed fixed step steepest descent method, which is repeated until 
no further progress is made; Step 3 is to relax the best solution gotten in Step 2 by a 
continuous optimization method conjugate gradient method. 

Doye et al ( |12J ) investigated the structures of clusters by mapping the structure of 
the global minimum as a function of both cluster size and the range of the pair potential 
which is appropriate to the clusters of diatomic molecule, Cgo molecule, and the ones 
between them both. For the larger clusters the structure of the global minimum changes 
from icosahedral to decahedral to fee as the range is decreased ( [12j ) . In ^3J Doye et al 
predicted the growth sequences for small decahedral and fee clusters by maximisation 
of the number of NN contacts. 

Lastly, in this section, we give a brief review of methods on the magic numbers 

= 17,23,24,72,88 and on the exceptions to [39]. Freeman et al ([E]) presented 
the best value for A^ = 17 when the thermodynamic properties of argon clusters were 
studied by a combination of classical and quantum Monte Carlo methods. The poly- 
icosahedral growth of Farges et al (|15j) starts from a 13-atom primitive icosahedron 
containing a central atom and 12 surface atoms. On each one of the five tetrahedral 
sites, surrounding a particular vertex, a new atom is added and finally a sixth atom 
is placed on top to create a pentagonal cap. In this way a 19-atom structure being 
made of double interpenetrating icosahedra, which is a 13-atom icosahedra sharing 9 
atoms, is obtained. I.e., for three pentagonal bipyramids each one shares an apex with 
its nearest neighbour. In this way a 23-atom model consisting of three interpenetrat- 
ing icosahedra is gotten for the best value known. Wille ([43]) used the simulated 
annealing method yielding low-lying energy states whose distribution depends on the 
cooling rate to find the best solution known for A^ = 24. Coleman et al ([9]) proposed 
a build-up process to construct the optimal solution structures. The HOC (half oc- 
tahedral cap) structure of the optimal solution for A^ = 72 is found by a prototype 
algorithm designed using the anisotropic effective energy simualted annealing method 
at each build-up stage ([9]). Wales and Doye ([33]) in 1997 gave the lowest values 
known for A^ = 192,201. Their method is so-called basin- hopping method, in which 
first the transformed function f{x) = min{/(a;)} was defined and performed starting 
from X by the PR conjugate gradient method (see, for example, (32]), and then the 
energy landscape for the function /(x) was explored using a canonical Monte Carlo 
simulation. Leary ([31]) has developed techniques for moving along sequences of local 
minima with decreasing energies to arrive at good candidates for global optima and got 
the best value known on A^ = 185. 
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3 Our approach to the Lennard-Jones Problem 



In this section we first briefly introduce the discrete gradient method; then present 
a hybrid discrete gradient and simulated annealing method. And then we give some 
techniques, with the discrete gradient method and the hybrid method, to tackle the 
Lennard-Jones Problem. 

3.1 Discrete Gradient Method 

Discrete gradient method is a derivative-free local search method for nonsmooth op- 
timization ([31 U]). Results of numerical experiments presented in [5] show that this 
method can jump over stationary points, which are not local minima, so we can reduce 
the number of stationary points, which we meet. Details of this method can be found, 
for example, in [3lll]. 

3.2 Hybrid Discrete Gradient and Simulated Annealing Method 

Simulated annealing method is a global search method. It has received a great deal 
of attention in last several years. First, this method was applied to combinatorial 
optimization, i.e. when the objective function is defined in a discrete domain (see 
[HI [30]). Later this method was applied to solve continuous global optimization problems 
(see [3 dni EHl [Ml [38]). Convergence of simulated annealing method is studied, for 
example, in [33j . 

Numerical experiments show that algorithms based on a combination of the global 
and local search techniques are more effective than the global search methods working 
alone (see, for example, [5l I48[ [20]). In these methods a local search algorithm is used 
to find a stationary point and a global search algorithm is used to escape from this 
stationary point and to find a new starting point for a local search algorithm. In this 
subsection, we develop a new hybrid simulated annealing and discrete gradient method. 
During the whole process of the method, both simulated annealing and discrete gradient 
methods need the the objective function values only. 

The simulated annealing method for solving continuous global optimization prob- 
lems was studied by many authors (see, for example, [3 [IHl [SHI [331 [38] ) ■ Neighborhood 
solution search procedure is one of the important steps in this method. Noticing that 
the neighborhood solution search for simulated annealing method should be at least 
based on two basic ideas: (a) neighbor means "nearby", (b) simulated annealing method 
is a stochastic method so that the neighborhood solution should be randomly taken, 
we may simply give a neighborhood solution search procedure for simulated annealing 
algorithm: 

Uniformly randomly keeping n — 1 coordinates of x, and making the left one coordinate of x 
randomly take a value such that the new solution x' still feasible. This gives x' . 

When the feasible region of the optimization problem is the unit simplex S = {x € 
M" : Y17=i Xi = I}, the neighborhood solution search procedure should be done a little 
modification: 

Uniformly randomly keeping n — 2 coordinates of x, and making one coordinate from the two 
elements left to x randomly take a value from [0,1] such that the value of the sum of the 
n — 1 coordinates gotten not greater than 1. Another left coordinate of x' is given the value 
1-sum. This gives x' . 
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The hybrid method starts from an initial point, first executes discrete gradient 
method to find local minimum, then carries on simulated annealing method in order 
to escape from this local minimum and to find a new starting point for the discrete 
gradient method. Then we again apply the discrete gradient method starting from the 
last point and so on until the sequence of the optimal objective function values gotten 
is convergent. The pseudo-code of the hybrid method is listed as following: 

Algorithm 1: Hybrid simulated annealing and discrete gradient method 
Initialization: 

Define the objective function f and its feasible solution space. 
Call the initial feasible solution generating procedure to get x. 
Call initial temperature selecting procedure to get T. 
Inialization of /: / = f(x). 

Initialize the neigbourhood feasible solution: xjneighbour = 0. 
Initialize xjyest: xjjest = x. 
Initialize f-best: fJbest = f. 

do { 

Discrete Gradient local search part: 

f -best-local = loca.l.sea.rch{x-best,X-new-gotten); 

X = xjnew_gotten\ 
Simulated Annealing global search part: 
do { 

do { 

xjneighbour = randomly 4)erturb{x); 

f -neighbour = f (xjneighbour); 

Calculate the difference A = fjneighbour — /; 

If (A < 0) or (random[0,l] < exp(-A/T)) 

X = xjneighbour f = fjneighbour; 
If (/ < f-best) xJbest = x f-best = f; 
} while (equilibrium has not been reached); 
Temperature annealing 
} while (Simulated Annealing stop criterion has not been met); 

} while ( f -best - f -best-local < -0.001 ); 

The convergence of the proposed hybrid method directly follows from the conver- 
gence of the simulated annealing and the discrete gradient methods. 

The simulated annealing method part consists of two procedures: inner procedure 
for the search of neighborhood solution and outer procedure for the cooling tempcratTire 
T. The number of iterations of inner and outer loops are taken to be large enough for 
guaranteeing that the sufficient iterations will be carried on to escape from the current 
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local minimum. In implementing the hybrid method, we use T = 0.9 *T as the temper- 
ature annealing schedule and the initial temperature is taken large enough according to 
the rule in [30]. We restrict the number of iterations for the outer procedure by 100 and 
number of iterations for the inner procedure by 1000. The discrete gradient method 
part is terminated when the distance between the approximation to the subdifferential 
and origin is less than a given tolerance e > (e = 10~^). 

3.3 Minimization of Lennard- Jones potential function 

In this subsection, the proposed discrete gradient method and hybrid discrete gradient 
and simulated annealing method, with the techniques given in the below, have been 
applied for minimization of well-known Lennard-Jones potential function in protein 
folding problem ([T]). The nonconvex Lennard-Jones potential function has a huge num- 
ber of local minima. Therefore many global optimization techniques fail to minimize 
it. The proposed hybrid method fails to solve this problem when number of atoms 
> 20. In order to reduce the number of local minima we suggest to approximate the 
function ^ ^ 

by the following function: 

c/(t) = max(5i(T),min(52(r),5'3(T))) 

where (71 (r) is the piecewise linear approximation of the function (/^(t) in segment (0, rg], 
52(7") is the piecewise linear approximation of this function over segment [ro,ri], and 
finally 53 (r) is the piecewise linear approximation over [ri , b] and b is large enough 
number. Here 



Such a approximation of the function y'(r) allows us to remove many local minima 
of the Lennard-Jones potential function and to get a good approximation to the global 
minimum of the objective function / in problem ([1]). 

In numerical experiments we take b = 16 and divide the segment [0.001, rg] into 100 
segments, the segment [ro, ri] into 100 segments and the [ri, 16] into 50 segments which 
allows one to get good approximations for the function </?(t). The replacement of the 
function ^{t) by the function g{T) makes the objective function nonsmooth. On the 
other side such a replacement significantly reduce the number of local minima. Since 
the discrete gradient method is a method of nonsmooth optimization the proposed 
hybrid method can be applied for solving this transformed problem. 

When solving the Lennard-Jones problem, first we use the discrete gradient method 
with build-up technique to relax to an initial solution. Then we apply the hybrid 
discrete gradient and simulated annealing method, with the good approximation for 
the objective function, to get another initial solution. Starting from this initial solution 
we again apply the derivative-free discrete gradient method and at last get the global 
solution. Results of numerical experiments are presented in Table [1] as the appendix. 
Results from Table [T] show that our techniques and methods effectively solve protein 
folding problem when number of atoms is not greater than 310. 

The structure of the optimal solutions corresponding to the above optimal values 
is illustrated and can be compared with the ones of [HI [39] in the following figures: 
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N=12, 13, 14, 15, 16, 17, 18, 19, 20 (ours) 
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N=12, 13, 14, 15, 16, 17, 18, 19, 20 ([49ii39j) 




21, 25, 27, 30, 34, 44 (ours) 




N=21, 25, 27, 30, 34, 44 ((491139 
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N=49, 56, 65, 67, 84, 93, 148, 170, 172 (ours) 




N=49, 56, 65, 67, 84, 93, 148, 170, 172 ([49,^) 



N=268, 288, 293, 298, 300, 301, 304, 308 (ours) 
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Table 1: Our numerical results for Lennard- Jones Protein Problem 
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-79 659789 
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-77 1 77043 


-77 1 77043 
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-81 684571 
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-86 573675 


-86 809789 


93 


-Q9 844461 


-Q9 844479 
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-97 34881 5 


-97 34881 5 






1 09 '^798*^*^ 
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-1 1 9 H9^^1 7 
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-1 98 096960 


-1 98 986571 


34 


-1 50 044598 


-1 50 044598 


44 


-907 631 655 


-907 688798 


4Q 


-939 091 863 


-93Q 091 864 




-983 '^94Q4'i 


-983 6431 05 
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-334 014007 


-334 971 539 


fi7 


-347 053308 


-347 959007 


84 


-459 96791 


-459 6573 


93 


-510 653123 


-510 8779 


148 


-881.072948 


-881.072971 


170 


-1024.791771 


-1024.791797 


172 


-1039.154878 


-1039.154907 


268 


-1706.182547 


-1706.182605 


288 


-1850.010789 


-1850.010842 


293 


-1888.427022 


-1888.427400 


298 


-1927.638727 


-1927.638785 


300 


-1942.106181 


-1942.106775 


301 


-1949.340973 


-1949.341015 


304 


-1971.044089 


-1971.044144 


308 


-1999.983235 


-1999.983300 



4 Conclusions 

In this paper, a brief review and a bibliography of important computational algorithms 
on minimizing the LJ potential energy are introduced in Sections 1 and 2. Section 3 
of this paper illuminates many beautiful graphs for the three dimensional structures of 
molecules with minimal LJ potential. 
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N=268, 288, 293, 298, 300, 301, 304, 308 (tMlES]) 
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